home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / SMOOTH11.ARJ / SMOOTH.C < prev    next >
C/C++ Source or Header  |  1991-06-10  |  22KB  |  793 lines
  1. /*    
  2.     NAME
  3.         smooth - split linear smoothing
  4.  
  5.     SYNOPSIS
  6.         smooth [file] [options]
  7.  
  8.         options are:
  9.           -a  [step [start]] automatic abscissas 
  10.           -b     break smooth after each label 
  11.           -c     general curve
  12.           -f  <num>  fraction of points to use for each fitted value 
  13.                  (default .5)
  14.           -n  <num>  number of points to use for each fitted value 
  15.                  (default 50%)
  16.           -r     print residuals rather than smoothed values
  17.           -s     split linear fit rather than "lowness"
  18.           -xl    take logs of x values before smoothing 
  19.           -yl    take logs of y values before smoothing 
  20.           -zl    take logs of z values before interpolating 
  21.                  (implies -3)
  22.           -3     3D case: x, y, and z given for each point
  23.  
  24.         Default smoothing method is "lowness": robust locally weighted 
  25.         regression & smoothing 
  26.  
  27.     METHODS
  28.         Lowness by W. S. Cleveland, and split linear fit by A. Owen
  29.  
  30.     AUTHOR
  31.         James Van Zandt (jrv@mbunix.mitre.org)
  32.  
  33. */
  34.  
  35. #include <stdio.h>
  36. #include <math.h>
  37.  
  38. #define VERSION "1.10"
  39.  
  40. #define MAX_WINDOWS 6
  41. #define MAX_POINTS 800
  42. #define MAXLABELS 25
  43.  
  44. int resid=0;                /* nonzero if residuals are requested */
  45. int fraction_specified=1;    /* nonzero if fit is to a specified fraction of 
  46.                                 the given data points */
  47. int number_specified=0;        /* nonzero if fit is to a specified number of 
  48.                                 the given data points */
  49. double fraction_fitted=.5;    /* fraction of points used for each fitted value */
  50. int num_fitted;                /* number of points used for each fitted value */
  51.  
  52. char *gmem();
  53.  
  54. /*
  55.     Split Linear Smoothing Algorithm
  56.  
  57.     Given:
  58.         A list of window sizes, SizeList, and n pairs (x[i],y[i]) sorted on x,
  59.  
  60.     Return:
  61.         the split linear smooth of y on x.
  62.  
  63.  
  64.         The general technique is due to Art Owen, who offers this discussion:
  65.  
  66.             "You should feel free to experiment with the algorithm,
  67.             since it has some ad hoc parts.  The essentials are: to use
  68.             uncentered windows of varying sizes along with the central
  69.             ones, to get zero weight on the worst fitting lines, and to
  70.             make the weight attached to a particular line size and
  71.             orientation vary smoothly as one traverses the data.  We
  72.             tried to find a simple way to meet all of these goals; the
  73.             algorithm we settled on was the simplest that worked for
  74.             us.
  75.                                 ...
  76.  
  77.             "West and Chan et.  al.  are useful for getting numerically
  78.             stable updating formulae for the regressions."
  79.  
  80.         references...
  81.  
  82.             John Alan McDonald and Art B. Owen, "Smoothing with Split
  83.             Linear Fits", LCS Technical Report No. 7, SLAC-PUB 3423,
  84.             AD-A149032, Laboratory for Computational Statistics, Dept. 
  85.             of Statistics, Stanford University, July 1984.
  86.  
  87.  
  88.             West, D.H.D., 1979, Updating Mean and Variance Estimates:
  89.             An Improved Method, Communications of the ACM, v 22, no.  9
  90.             p 532-535 (1979).
  91.  
  92.             Chan, T.F., Golub, G.H., and Leveque, R.J., 1983,
  93.             Algorithms for Computing the Sample Variance: Analysis and
  94.             Recommendations, The American Statistician v 37, p 242-247
  95.             (1983).
  96.  
  97.     IMPLEMENTATION
  98.         For general curves, the given x- and y- (and z-, if present)
  99.         points are regarded as functions of the distance along a
  100.         smoothed path.
  101.  
  102.         The arrays take a lot of space.  For n points, the number of
  103.         doubles is approximately 38*n, plus 2*n for general curve, plus
  104.         2*n for 3D case.  For 100 points and 8 byte doubles, this means
  105.         at least 8*38*100=30400 bytes.
  106.  
  107.         Execution time...  101 points took 151 seconds on 7.5 MHz V-20,
  108.         no 8087.
  109.  
  110.  
  111.         The updating formulas mentioned above are not used in this
  112.         program.  The selection of window sizes (a geometric sequence)
  113.         is my own.  -JVZ
  114.         
  115. */
  116.  
  117. SplitLinearSmooth(SizeList,x,y,Smooth) int SizeList; double Smooth[],x[],y[];
  118. {    int i;
  119.     BasicSplitLinearSmooth( SizeList, x, y, Smooth );
  120.     BasicSplitLinearSmooth( SizeList, x, Smooth, Smooth );
  121.     if(resid)
  122.         for (i=0; i<SizeList; i++)
  123.             Smooth[i] = y[i]-Smooth[i];
  124. }
  125.  
  126. /*
  127.     Basic Split Linear Smooth  Algorithm    by Art Owen  <owen@su-score>
  128.  
  129.     Given:
  130.         A list of window sizes and n pairs (x[i],y[i]) sorted on x,
  131.  
  132.     Return:
  133.         the basic split linear smooth of y on x.
  134. */
  135. #define f(x,y,z) ff[((x)*sizes+(y))*3+(z)]
  136. #define pmse(x,y,z) errors[((x)*sizes+(y))*3+(z)]
  137. #define left    0
  138. #define center    1
  139. #define right    2
  140.  
  141. static double *ff, *errors;
  142.  
  143. BasicSplitLinearSmooth(n, x, y, Smooth) int n; double x[], y[], Smooth[];
  144. {    int k, C, L, R, i, sizes, num, start, finish, o, ws;
  145.     int window_size[MAX_WINDOWS]; 
  146.     double err, del, alpha, beta, w_total, w, s, c;
  147.     double sum, sumx, sumxx, sumy, sumxy;
  148.  
  149.         /*
  150.             window sizes must be _odd_ in this implementation
  151.             sizes are: 5, 11, 23, 47, ... 6*2^k-1
  152.             or else:   7, 15, 31, 63, ... 8*2^k-1 
  153.             depending on the maximum number of points to be fitted
  154.         */
  155.     for (i = num_fitted+1; i >= 8; i/=2) {}
  156.     if(i >= 6) i=6; else i=8;
  157.     for (sizes=0; i-1 <= num_fitted && sizes < MAX_WINDOWS; i *= 2)
  158.         window_size[sizes++]=i-1;
  159.  
  160.             /* most of the memory is allocated by these two statements */
  161.     ff=(double *)gmem(n*sizes*3*sizeof(double));
  162.     errors=(double *)gmem(n*sizes*3*sizeof(double));
  163.     for (k=0; k<sizes; k++)
  164.         {sum=sumx=sumxx=sumy=sumxy=0.;
  165.         /* Initialize the regression-window data structure */
  166.         ws=window_size[k];
  167.         C = -(ws/2);
  168.         L = C - ws/2;
  169.         R = C + ws/2;
  170.         for (  ; L<n ; L++, C++, R++) 
  171.             {
  172.             if ( R < n )
  173.                 {/* update the regression-window with (x[R],y[R]) */
  174.                 /* Update routines are described in West (1979) */
  175.                 /* these are my guesses -JVZ */
  176.                 sum+=1.; sumx+=x[R]; sumxx+=x[R]*x[R];
  177.                 sumy+=y[R]; sumxy+=x[R]*y[R];
  178.                 }
  179.             if ( L-1 >= 0 )
  180.                 {/* downdate the regression-window with (x[L-1],y[L-1]) */
  181.                 sum-=1.; sumx-=x[L-1]; sumxx-=x[L-1]*x[L-1];
  182.                 sumy-=y[L-1]; sumxy-=x[L-1]*y[L-1];
  183.                 }
  184.             if ( sum>=ws ) /* the window covers enough points */
  185.                 {/* Fit linear regression coefs. alpha, beta  */
  186.                 alpha=(sumy*sumxx - sumx*sumxy)/(sum*sumxx - sumx*sumx);
  187.                 beta=(sumxy - alpha*sumx)/sumxx;
  188.                 if ( 0 <= R && R < n )
  189.                     {f(R,k,left)= alpha + beta*x[R];
  190.                     /* compute pmse(R,k,left) */
  191. /*
  192. pmse(i,k,o) = 
  193.     ( sum{j in W(i,k,o)-{i}} (y[j]-alpha-beta*x[j])^2) / ( |W(i,k,o)|-3 )
  194. where W is the window of points.
  195. */
  196.                     start=R-ws;
  197.                     if(start<0) start=0;
  198.                     if(R-start>3)
  199.                         {err=0.;
  200.                         for (i=start; i<R; i++)
  201.                             {del = y[i]-alpha-beta*x[i];
  202.                             err += del*del;
  203.                             }
  204.                         pmse(R,k,left) = err/(R-start-3);
  205.                         }
  206.                     }
  207.                  if ( 0 <= C && C < n )
  208.                     {f(C,k,center)= alpha + beta*x[C]; 
  209.                     /* compute pmse(C,k,center) */
  210.                     start=C-ws/2;
  211.                     if(start<0) start=0;
  212.                     finish=C+ws/2+1;
  213.                     if(finish>n) finish=n;
  214.                     if(finish-start>3)
  215.                         {err=0.;
  216.                         for (i=start; i<finish; i++)
  217.                             {if(i==C) continue;
  218.                             del = y[i]-alpha-beta*x[i];
  219.                             err += del*del;
  220.                             }
  221.                         pmse(C,k,center) = err/(finish-start-3);
  222.                         }
  223.                     }
  224.                 if ( 0 <= L )
  225.                     {f(L,k,right)= alpha + beta*x[L];
  226.                     /* compute pmse(L,k,right) */
  227.                     err=0.;
  228.                     finish=L+ws;
  229.                     if(finish>n) finish=n;
  230.                     if(finish-L>3)
  231.                         {for (i=L; i<finish; i++)
  232.                             {del = y[i]-alpha-beta*x[i];
  233.                             err += del*del;
  234.                             }
  235.                         pmse(L,k,right) = err/(finish-L-3);
  236.                         }
  237.                     }
  238.                 }
  239.             else    /* window doesn't cover enough points */
  240.                 {if( R<n ) pmse(R,k,left) = -1;
  241.                 if( 0<=C && C<n ) pmse(C,k,center) = -1;
  242.                 if( 0<=L ) pmse(L,k,right) = -1;
  243.                 }
  244.             }
  245.         }
  246.     for (i=0; i<n; i++)
  247.         {/* pmse cutoff: c = Average{k,o} pmse(i,k,o) */
  248.         err=0.;
  249.         num=0;
  250.         for (k=0; k<sizes; k++)
  251.             for (o=0; o<3; o++)
  252.                 if(pmse(i,k,o)>=0.) {err += pmse(i,k,o); num++;}
  253.         c=err/num;
  254.         s = w_total = 0.;
  255.         for (k=0; k<sizes; k++)
  256.             for (o=0; o<3; o++)
  257.                 {err=pmse(i,k,o);
  258.                 if(err>=0. && err<c) /* pmse is present, and below cutoff */
  259.                     {w = (err-c)*(err-c);
  260.                     s += w*f(i,k,o);
  261.                     w_total += w;
  262.                     }
  263.                 }
  264.             /* Smooth[i] = sum{k,o} w(i,k,o)*f(i,k,o) / sum{k,o} w(i,k,o) */
  265.         Smooth[i]=s/w_total;
  266.         }
  267.     free(errors);
  268.     free(ff);
  269. }
  270.  
  271. double abscissa=0.;
  272. double abscissa_step=1.;
  273. double *x, *y, *z, *xs, *ys, *zs, *path;
  274.  
  275. int xlog=0;        /* nonzero if taking log of x values */
  276. int ylog=0;        /* nonzero if taking log of y values */
  277. int zlog=0;        /* nonzero if taking log of z values */
  278. int debugging=0;
  279. int split=0;    /* nonzero if split linear smooth is requested */
  280. int breaking=0;        /* nonzero if breaking at labels */
  281. int labels=0;        /* number of labels in input data */
  282. int threed=0;    /* nonzero if 3D case (x, y, and z given for each point) */
  283. int automatic_abscissas=0;
  284. int abscissa_arguments=0;
  285. int curve=0;        /* nonzero if general curve permitted */
  286. int x_arguments=0;
  287. int n;                /* number of MAX_POINTS in x, y, and ys */
  288. int index_array[MAXLABELS];    /* indexes into x and y */
  289. int *p_data=index_array;
  290. int total=100;
  291.  
  292. FILE *ifile;
  293.  
  294. char *label_array[MAXLABELS];
  295. char **p_text=label_array;
  296.  
  297. main(argc,argv) int argc; char **argv;
  298. {    int nn,origin;
  299.     read_data(argc,argv);
  300.  
  301.                             /* check switch arguments */
  302.     if(fraction_specified)
  303.         {if(fraction_fitted<0. || fraction_fitted>1.)
  304.             {fprintf(stderr, "fraction fitted=%f must be in range 0 to 1\n", 
  305.                                                             fraction_fitted);
  306.             exit(1);
  307.             }
  308.         num_fitted = floor(fraction_fitted*n+.5);
  309.         }
  310.     nn = 1;
  311.     if(split) nn = 11;
  312.     if(num_fitted < nn || num_fitted > n)
  313.         {fprintf(stderr, "number of points fitted=%d must be in range %d...%d", 
  314.                                                         num_fitted, nn, n);
  315.         exit(1);
  316.         }
  317.  
  318.     if(breaking && labels)
  319.         {origin=0;
  320.         while(labels--)
  321.             {p_data[0] -= origin;
  322.             nn=p_data[0]+1;
  323.             if(nn) curv0(nn,total);
  324.             origin += nn;
  325.             path += nn;
  326.             x += nn;
  327.             y += nn;
  328.             p_data++;
  329.             p_text++;
  330.             }
  331.         }
  332.     else
  333.         {curv0(n,total);
  334.         }
  335.     exit(0);
  336. }
  337.  
  338. curv0(n,requested)
  339. int n,requested;
  340. {    int i, j, each, stop, seg=0;
  341.     double s, ds, xx, yy, zz;
  342.     if(n<5)                        /* just copy a short file */
  343.         {for (j=0; j<n; j++)
  344.             {if(!automatic_abscissas)
  345.                 {xx=x[n-1]; if(xlog) xx=exp(xx);
  346.                 printf("%g ", xx);
  347.                 }
  348.             yy=y[n-1]; if(ylog) yy=exp(yy);
  349.             printf("%g ", yy);
  350.             if(threed)
  351.                 {zz=z[n-1]; if(zlog) zz=exp(zz);
  352.                 printf("%g ",zz);
  353.                 }
  354.             if(j==p_data[seg]) printf(p_text[seg++]);
  355.             putchar('\n');
  356.             }
  357.         return;
  358.         }
  359.     if(split)
  360.         {if(curve)
  361.             {SplitLinearSmooth( n, path, x, xs);
  362.             SplitLinearSmooth( n, path, y, ys);
  363.             if(threed) SplitLinearSmooth(n, path, z, zs);
  364.             }
  365.         else
  366.             {SplitLinearSmooth( n, x, y, ys);
  367.             if(threed) SplitLinearSmooth( n, x, z, zs);
  368.             }
  369.         }
  370.     else /* lowness smooth */
  371.         {if(curve)
  372.             {lowness( n, path, x, xs);
  373.             lowness( n, path, y, ys);
  374.             if(threed) lowness(n, path, z, zs);
  375.             }
  376.         else
  377.             {lowness( n, x, y, ys);
  378.             if(threed) lowness( n, x, z, zs);
  379.             }
  380.         }
  381.     for (j=0; j<n; j++)
  382.         {if(!automatic_abscissas)
  383.             {xx = curve ? xs[j] : x[j];
  384.             printf("%g ", xlog ? exp(xx) : xx);
  385.             }
  386.         printf("%g ", ylog ? exp(ys[j]) : ys[j] );
  387.         if(threed) printf("%g ", zlog ? exp(zs[j]) : zs[j] );
  388.         if(j==p_data[seg]) printf(p_text[seg++]);
  389.         putchar('\n');
  390.         }
  391. }
  392.  
  393. double mylog(x) double x;
  394. {    if(x<=0.)
  395.         {fprintf(stderr, "can\'t take log of %f\n", x);
  396.         exit(1);
  397.         }
  398.     return log(x);
  399. }
  400.  
  401. read_data(argc,argv) int argc; char **argv;
  402. {    int i, j, length, skipping, ac;
  403.     double xx, yy, zz, d, *pd, sum, xp, yp, zp, xq, yq, zq;
  404.     char *s,*t, **av, *strchr();
  405. #define BUFSIZE 256
  406.     static char buf[BUFSIZE+2];
  407.  
  408.     argc--; argv++;
  409.     ac=argc; av=argv; argc=0;
  410.     while(ac>0)
  411.         {if(**av=='?') help();
  412.         else if(**av=='-' && (i=get_parameter( ac, av )) && i > 0)
  413.             {ac-=i; av+=i;
  414.             }
  415.         else {argv[argc++] = *av++; ac--;}
  416.         }
  417.  
  418.     if(argc<1 || **argv=='-') ifile=stdin;
  419.     else
  420.         {if(argc>=1)
  421.             {ifile=fopen(argv[0],"r");
  422.             if(ifile==NULL) {printf("file %s not found\n",argv[0]); exit(1);}
  423.             argc--; argv++;
  424.             }
  425.         else help();
  426.         }
  427.     x=path=(double *)gmem(MAX_POINTS*sizeof(double));
  428.     y=(double *)gmem(MAX_POINTS*sizeof(double));
  429.     if(threed) z=(double *)gmem(MAX_POINTS*sizeof(double));
  430.     p_data[0]=-1;
  431.     skipping=2;
  432.     if(automatic_abscissas) skipping--;
  433.     if(threed) skipping++;
  434.     i=0;
  435.     while(i<MAX_POINTS)
  436.         {if(fgets(buf,BUFSIZE,ifile)==0)break;
  437.         buf[strlen(buf)-1]=0;                    /* zero out the line feed */
  438.         if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
  439.         if(t = strchr(buf, ';')) *t=0;            /* zap same-line comment */
  440.         for (t=buf; *t && isspace(*t); t++) {}    /* find first nonblank char */
  441.         if(*t == 0) continue;                    /* skip blank lines */
  442.         if(automatic_abscissas)
  443.             {x[i]=abscissa;
  444.             abscissa+=abscissa_step;
  445.             if(threed) sscanf(buf,"%lf %lf",&y[i],&z[i]);
  446.             else sscanf(buf,"%lf",&y[i]);
  447.             if(ylog) y[i]=mylog(y[i]);
  448.             if(zlog) z[i]=mylog(z[i]);
  449.             }
  450.         else
  451.             {if(threed) sscanf(buf,"%lf %lf %lf",&x[i],&y[i],&z[i]);
  452.             else sscanf(buf,"%lf %lf",&x[i],&y[i]);
  453.             if(xlog) x[i]=mylog(x[i]);
  454.             if(ylog) y[i]=mylog(y[i]);
  455.             if(zlog) z[i]=mylog(z[i]);
  456.             }
  457.         s=buf;                      /* start looking for label */
  458.         for (j=skipping; j--; )
  459.             {while(*s==' ')s++;            /* skip a number */
  460.             while(*s && (*s!=' '))s++;
  461.             }
  462.         while(*s==' ')s++;
  463.         if((length=strlen(s))&&(labels<MAXLABELS))
  464.             {if(*s=='\"')           /* label in quotes */
  465.                 {t=s+1;
  466.                 while(*t && (*t!='\"')) t++;
  467.                 t++;
  468.                 }
  469.             else                    /* label without quotes */
  470.                 {t=s;
  471.                 while(*t && (*t!=' '))t++;
  472.                 }
  473.             *t=0;
  474.             length=t-s;
  475.             p_data[labels]=i;
  476.             p_text[labels]=gmem(length+1);
  477.             if(p_text[labels]) strcpy(p_text[labels++],s);
  478.             }
  479.         i++;
  480.         }
  481.     n=i;
  482.     if(n <= 1) {fprintf(stderr, "too little data"); exit(1);}
  483.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  484.         {p_data[labels]=i-1;
  485.         if(p_text[labels]=gmem(1)) *p_text[labels++]=0;
  486.         }
  487.     ys=(double *)gmem(n*sizeof(double));
  488.     if(threed) zs=(double *)gmem(n*sizeof(double));
  489.     if(curve)
  490.         {xs=(double *)gmem(n*sizeof(double));
  491.         path=(double *)gmem(n*sizeof(double));
  492.         path[0]=sum=0.;
  493.         if(threed)
  494.             {xp=4.*x[0]; yp=4.*y[0]; zp=4.*z[0];
  495.             for (i=1; i<n; i++)
  496.                 {xq=x[i-1]+2.*x[i]+x[i+1];
  497.                 yq=y[i-1]+2.*y[i]+y[i+1];
  498.                 zq=z[i-1]+2.*z[i]+z[i+1];
  499.                 xx=xq-xp;
  500.                 yy=yq-yp;
  501.                 zz=zq-zp;
  502.                 path[i]=(sum+=sqrt(xx*xx + yy*yy + zz*zz));
  503.                 xp=xq; yp=yq; zp=zq;
  504.                 }
  505.             xq=4.*x[n-1]; yq=4.*y[n-1]; zq=4.*y[n-1];
  506.             xx=xq-xp;
  507.             yy=yq-yp;
  508.             zz=zq-zp;
  509.             path[i]=(sum+=sqrt(xx*xx + yy*yy + zz*zz));
  510.             }
  511.         else
  512.             {xp=4.*x[0]; yp=4.*y[0];
  513.             for (i=1; i<n-1; i++)
  514.                 {xq=x[i-1]+2.*x[i]+x[i+1];
  515.                 yq=y[i-1]+2.*y[i]+y[i+1];
  516.                 xx=xq-xp;
  517.                 yy=yq-yp;
  518.                 path[i]=(sum+=sqrt(xx*xx + yy*yy));
  519.                 xp=xq; yp=yq;
  520.                 }
  521.             xq=4.*x[n-1]; yq=4.*y[n-1];
  522.             xx=xq-xp;
  523.             yy=yq-yp;
  524.             path[i]=(sum+=sqrt(xx*xx + yy*yy));
  525.             }
  526.         }
  527.     else if(split)
  528.         {for(i=1; i<n; i++)
  529.             {if(x[i]<=x[i-1])
  530.                 {fprintf(stderr,
  531.                 "ERROR: independent variable not strictly increasing\n");
  532.                 exit(1);
  533.                 }
  534.             }
  535.         }
  536.     else
  537.         {for(i=1; i<n; i++)
  538.             {if(x[i]<x[i-1])
  539.                 {fprintf(stderr,
  540.                 "ERROR: data must be sorted on x value\n");
  541.                 exit(1);
  542.                 }
  543.             }
  544.         }
  545. }
  546.  
  547.  
  548. /* get_parameter - process one command line option
  549.         (return # parameters used) */
  550. get_parameter(argc,argv) int argc; char **argv;
  551. {    int i, permitted;
  552.     double temp;
  553.     if(strcmp(*argv,"-a")==0)
  554.         {automatic_abscissas=1;
  555.         if(argc<2 || *argv[1]=='-') return 1; /* no negative step */
  556.         permitted=2;
  557.         if(argc<3 || *argv[2]=='-') permitted=1; /* no neg. starting abscissa */
  558.         i=get_double(argc,argv,permitted,&abscissa_step,&abscissa,&abscissa);
  559.         abscissa_arguments=i-1;
  560.         return i;
  561.         }
  562.     else if(strcmp(*argv, "-f")==0)
  563.         {i=get_double(argc, argv, 1, &fraction_fitted, &fraction_fitted, &fraction_fitted);
  564.         fraction_specified=1;
  565.         number_specified=0;
  566.         return i;
  567.         }
  568.     else if(strcmp(*argv, "-n")==0)
  569.         {i=get_double(argc, argv, 1, &temp, &temp, &temp);
  570.         fraction_specified=0;
  571.         number_specified=1;
  572.         num_fitted = temp;
  573.         return i;
  574.         }
  575.     else if(strcmp(*argv, "-b")==0) {breaking=1;        return 1;}
  576.     else if(strcmp(*argv, "-c")==0) {curve=1;           return 1;}
  577.     else if(strcmp(*argv, "-d")==0) {debugging=1;       return 1;}
  578.     else if(strcmp(*argv, "-r")==0) {resid=1;           return 1;}
  579.     else if(strcmp(*argv, "-s")==0) {split=1;           return 1;}
  580.     else if(strcmp(*argv,"-xl")==0) {xlog++;            return 1;}
  581.     else if(strcmp(*argv,"-yl")==0) {ylog++;            return 1;}
  582.     else if(strcmp(*argv,"-zl")==0) {zlog++; threed++;  return 1;}
  583.     else if(strcmp(*argv, "-3")==0) {threed++;          return 1;}
  584.     else gripe(argv);
  585. }
  586.  
  587. get_double(argc,argv,permitted,a,b,c)
  588. int argc,permitted; char **argv; double *a,*b,*c;
  589. {    int i=1;
  590.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  591.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  592.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  593.     return i;
  594. }
  595.  
  596. gripe_arg(s) char *s;
  597. {    fprintf(stderr,"argument missing for switch %s",s);
  598.     help();
  599. }
  600.  
  601. gripe(argv) char **argv;
  602. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  603.     help();
  604. }
  605.  
  606. char *msg[]=
  607. {"smooth   version ", VERSION, "\n",
  608. "\nusage: smooth [file] [options]\n",
  609. "options are:\n",
  610. "  -a  [step [start]] automatic abscissas \n",
  611. "  -b     break smooth after each label \n",
  612. "  -c     general curve \n",
  613. "  -f  <num>  fraction of points to use for each fitted value (default .5)\n",
  614. "  -n  <num>  number of points to use for each fitted value (default 50%%)\n",
  615. "  -r     print residuals rather than smoothed values\n",
  616. "  -s     split linear fit rather than \"lowness\"\n",
  617. "  -xl    take logs of x values before smoothing \n",
  618. "  -yl    take logs of y values before smoothing \n",
  619. "  -zl    take logs of z values before interpolating \n",
  620. "         (implies -3)\n",
  621. "  -3     3D case: x, y, and z given for each point\n",
  622. "\n",
  623. "Default smoothing method is \"lowness\": robust locally weighted \n",
  624. "regression & smoothing \n",
  625. 0
  626. };
  627.  
  628. help()
  629. {    char **s;
  630.     s=msg;
  631.     while(*s) printf(*s++);
  632.     exit(1);
  633. }
  634.  
  635. numeric(s) char *s;
  636. {    char c;
  637.     while(c=*s++)
  638.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  639.         return 0;
  640.         }
  641.     return 1;
  642. }
  643.  
  644. char *gmem(size) unsigned int size;
  645. {    char *p;
  646.     void *malloc();
  647.     p=(char *)malloc(size);
  648.     if(!p)
  649.         {fprintf(stderr,"can\'t allocate temporary storage of %u\n",size);
  650.         exit(1);
  651.         }
  652.     return p;
  653. }
  654.  
  655. fpcompare(r, s) double *r, *s;
  656. {    double d;
  657.     d=*r-*s;
  658.     if(d<0.) return -1;
  659.     if(d>0.) return 1;
  660.     return 0;
  661. }
  662.  
  663. /*
  664.     Lowness
  665.  
  666.         Robust locally weighted regression is a method for smoothing a
  667.         scatterplot, (x[i], y[i]), i=1,...,n, in which the fitted value
  668.         at x[k] is the value of a polynomial fit to the data using
  669.         weighted least squares, where the weight for (x[i], y[i]) is
  670.         large if x[i] is close to x[k].  Robustness is added by
  671.         calculating residuals and repeating the procedure with reduced
  672.         weights on points with large residuals.
  673.  
  674.     Reference: 
  675.         W.  S.  Cleveland, "Robust Locally Weighted Regression and
  676.         Smoothing Scatterplots", Journal of the American Statistical
  677.         Association, v74, n368, p829 (Dec 79)
  678.  
  679. */
  680.  
  681. lowness(n, x, y, ys) int n; double *x, *y, *ys;
  682. {    double *t, *r, *w, d, d1, u, m, m6i, sum, a[2];
  683.     int i, j, k;
  684.  
  685.     num_fitted=floor(fraction_fitted*n+.5);
  686.     t=(double *)gmem(num_fitted*sizeof(double));
  687.     r=w=(double *)gmem(n*sizeof(double)); /* r and w share space */
  688.     k=0;
  689.     for (i=0; i<n; i++)
  690.         {if(i==0 || x[i]!=x[i-1])
  691.             {while(k+num_fitted<n && (x[i]-x[k] > x[k+num_fitted]-x[i])) k++;
  692.             d=x[i]-x[k]; d1=x[k+num_fitted-1]-x[i];
  693.             if(d1>d) d=d1;
  694.             if(d==0.)
  695.                 {sum=0.; 
  696.                 for (j=0; j<num_fitted; j++) sum += y[i+j];
  697.                 a[0]=sum/num_fitted; a[1]=0.;
  698.                 }
  699.             else
  700.                 {for (j=0; j<num_fitted; j++)     /* calculate tricube weighting */
  701.                     {u=fabs(x[i]-x[k+j])/d; u=1.-u*u*u; t[j]=u*u*u;
  702.                     }
  703.                 lin(num_fitted, x+k, t, y+k, a);
  704.                 }
  705.             }
  706.         ys[i]=r[i]=fabs(y[i]-(a[0]+a[1]*x[i]));
  707.         }
  708.     /* sorting to find median value - although it's overkill */
  709.     qsort(ys, n, sizeof(double), fpcompare);
  710.     m=ys[n/2];
  711.     if(m==0.)    /* input function is already linear - can't smooth any more */
  712.         {for (i=0; i<n; i++) ys[i]=y[i];
  713.         free(t); free(w);
  714.         return;
  715.         }
  716.     m6i=1./(6.*m);
  717.     for (i=0; i<n; i++)
  718.         {u=r[i]*m6i;
  719.         if(u<1.) {u=1.-u*u; w[i]=u*u;}
  720.         else w[i]=0.;
  721.         }
  722.     k=0;
  723.     for (i=0; i<n; i++)
  724.         {if(i==0 || x[i]!=x[i-1])
  725.             {while(k+num_fitted<n && (x[i]-x[k] > x[k+num_fitted]-x[i])) k++;
  726.             d=x[i]-x[k]; d1=x[k+num_fitted-1]-x[i];
  727.             if(d1>d) d=d1;
  728.             if(d==0.)
  729.                 {sum=0.; 
  730.                 for (j=0; j<num_fitted; j++) sum += y[i+j];
  731.                 a[0]=sum/num_fitted; a[1]=0.;
  732.                 }
  733.             else
  734.                 {for (j=0; j<num_fitted; j++)     /* calculate tricube weighting */
  735.                     {u=fabs(x[i]-x[k+j])/d; u=1.-u*u*u; t[j]=u*u*u*w[k+j];
  736.                     }
  737.                 lin(num_fitted, x+k, t, y+k, a);
  738.                 }
  739.             }
  740.         ys[i]=a[0] + a[1]*x[i];
  741.         if(resid) ys[i]=y[i]-ys[i];
  742.         }
  743.     free(t); free(w);
  744. }
  745.  
  746. /*    lin - linear least squares fit */
  747.  
  748. lin(n, x, w, y, a) 
  749. int n;             /* # points in x, w, or y */
  750. double *x,         /* independent variables */
  751. *w,             /* weights */
  752. *y,             /* dependent variables */
  753. *a;                /* resulting coeficients... y[i] = a[1]*x[i] + a[0] (approx) */
  754. {    int i, j;
  755.     double sum, sumx, sumxx, sumy, sumxy, xw, cond, decomp();
  756.     static double d[2][2];
  757.     static int ipvt[2];
  758.  
  759.     sum=sumx=sumxx=sumy=sumxy=0.;
  760.     for (i=0; i<n; i++)
  761.         {sum += *w;
  762.         sumx += (xw = *x * *w);
  763.         sumxx += *x * xw;
  764.         sumy += *y * *w;
  765.         sumxy += *y * xw;
  766.         x++; y++; w++;
  767.         }
  768.     d[0][0]=sum; d[0][1]=sumx; a[0]=sumy;
  769.     d[1][0]=sumx;  d[1][1]=sumxx;  a[1]=sumxy;
  770.     cond=decomp(2, 2, d, ipvt);
  771.     if(cond==(cond+1.)) 
  772.         {fprintf(stderr, "ill conditioned matrix"); 
  773. /********
  774. **********/
  775. showv("x=", x-n, n);
  776. showv("w=", w-n, n);
  777. showv("y=", y-n, n);
  778. printf("cond=%9.3e\n", cond);
  779.         exit(1);
  780.         }
  781.     solve(2, 2, d, a, ipvt);
  782. }
  783.  
  784. showv(s,a,n) char *s; double a[]; int n;
  785. {    int i,j;
  786.     printf("%s \n",s);
  787.     for (i=0; i<n; i++)
  788.         {printf("%8.4f ",a[i]);
  789.         if(i%8==7) printf("\n");
  790.         }
  791.     printf("\n");
  792. }
  793.